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ABSTRACT 


The ability to determine the structural dynamics of space- 
based platforms from ground-based radar resolved Doppler 
measurements will aid in the study of control/structure 
interaction. The Naval Research Laboratory and Lincoln 
Laboratory conducted an experiment to determine the 
feasibility of this method. To accomplish this experiment the 
LACE satellite was equipped with retroreflectors and the 
ground-based Firepond laser radar facility was employed. 
Vibrational information is found from the difference between 
the reflected Doppler frequencies of the retroreflectors. The 
method of extracting the Doppler separation was to obtain the 
power spectrum of the heterodyne signal envelope. A pulse-by- 
pulse processing of the data yields the Doppler separation 
history over time. Due to a relatively large amount of 
clutter in the processed data, a filtering mechanism was 
employed. The histogram technique is the current filtering- 
based method employed to obtain a Doppler separation history. 
This thesis addresses the implementation of the Kalman filter 
algorithm in conjunction with the Rauch-Tung-Striebel fixed- 
interval optimal smoother algorithm to perform this filtering 
task. The Kalman smoother filtering based method of processing 
the data produced superior results when compared with the 


histogram filtering based method. 
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I. INTRODUCTION 


Control/structure interaction (CSI) is the interaction 
between control systems and the platform or the structural 
appendages. The problem of control/structure interaction is 
of great interest in the development of new space-based 
platforms. There are unique problems associated with the 
ground-based testing of structures designed for weightless 
environments. In strong gravitational fields spaced-based 
platforms exhibit different structural characteristics than 
those found in a weightless environment. This presents 
problems in developing models to simulate the structural 
dynamics based on data obtained from ground level 
experimentation. An ideal way to develop accurate models using 
experimental data is to obtain the data while the platform is 
in orbit. There are various methods that could be utilized to 
accomplish this task. A method that does not involve 
telemetric links or sophisticated electronic hardware 
installed on the platform is remote ground-based Doppler 
resolved measurements. These ground-based Doppler resolved 
measurements will then be used to analyze the structural 
dynamics of the platform. The Naval Research Laboratory [Ref. 
1] and Lincoln Laboratory [{Ref. 2] have been sponsoring a 
series of experiments to determine the feasibility of using 


the ground-based Doppler resolved measurement approach. 


To accomplish this study, the Laser Atmospheric 
Compensated Experiment (LACE) Satellite (object number 20496) 
was equipped with three germanium IR retroreflectors prior to 
launch. These IR retroreflectors were located on the forward 
boom, trailing boom, and body of the satellite as indicated in 


Figure 1. 
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Figure 1. LACE Spacecraft [Fig la from Ref. 1] 


The retroreflectors are then illuminated with the ground- 
based Firepond laser radar as depicted in Figure 2. The 


Firepond is a coherent narrowband 10.6 micrometer laser radar 
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Figure 2. Primary Experiment {Fig lb from Ref. 1] 


facility. The reflected signals from the satellite contained 


Doppler frequency components proportional to the relative 


velocity of the germanium retroreflectors projected along the 


radar line-of-sight. 
between the Doppler 
procedure employed 


separation from the 


The Doppler separation is the difference 
frequencies of the retroreflectors. The 
in the extraction of this Doppler 


reflected signal consisted of obtaining 


the power spectrum of the heterodyne signal envelope. This 


method greatly reduced the tolerance requirement of the 


equipment and calculations needed to extract the Doppler 


separation data from the measurements. 


One major drawback to this method is that noise rejection 
for this system is poor. The Doppler separation history 
obtained from the pulse-by-pulse processing had a significant 
amount of clutter. Lincoln Laboratory used a histogram 
filtering based technique for further noise rejection to 
obtain a Doppler separation history. Another approach to this 
problem, explored in the following chapters, will consist of 
the use of a Kalman Filter in conjunction with the 
Rauch-Tung-Striebel fixed interval optimal smoother. The 
Kalman filter is used to estimate the Doppler frequency and 
control a dynamic tracking window. This method of extracting 
the information components from the data showed a remarkable 
improvement in the resolution of the LACE satellite's Doppler 
separation history over the histogram filtering-based 
technique. 

The basic theory and the equations used in the LACE 
signal processing are discussed in the second chapter. In 
Chapter III the Kalman filter equations and performance 
characteristics are discussed. The Rauch-Tung-Striebel fixed 
interval smoother is discussed in Chapter IV. The remaining 
two chapters are devoted to the description of how the Kalman 
filter and the Rauch-Tung-Striebel fixed interval optimal 
smoother were utilized in this application and the results 


obtained. 


A. THE LACE DYNAMICS EXPERIMENT 

Lincoln Laboratory discusses the experiment and the 
procedures employed for the analysis of the narrowband IR 
measurements obtained for the 13 and 18 July 1990 [Ref. 2]. On 
both days, the LACE satellite's configuration had the leading 
boom's extension at 4.6 meters (15 feet) and the trailing 
boom's extension at 46 meters (150 feet). This configuration 
had been set up for a Significant time prior to the 
illumination to eliminate any vibrational modes that might 
have been excited by the boom's movement. On both tracking 
runs, the Firepond laser radar had a peak transmit power of 
780 watts, a pulse duration of 1.5 milliseconds, and a pulse 
repetition frequency (PRF) of 62.5 Hz. The maximum elevation 
that the LACE satellite achieved relative to the Firepond 
laser radar site on 13 July was 82 degrees and on 18 July was 
77 degrees. Due to the transmission beam's footprint of 12 
meters at the minimum range of 547 kilometers, only the 
leading boom's retroreflector and body's retroreflector were 
illuminated. 

The procedure that was used in the analysis of the 
received data first consisted of digitizing 3.4 milliseconds 
of the in-phase and quadrature (IQ) data at 1.2 MHz (generating 
4080 complex IQ samples). These complex IQ samples were then 
squared to yield the IQ envelope as shown in Figure 3. This 
figure shows the first radar return obtained from the tape for 


18 July 1990. The power spectral density from the IQ envelope 


was then obtained. Figure 4 depicts the power spectral 
density of Figure 3. Observation of this radar return 
indicates a nominal 12 kHz Doppler separation. 

The Doppler separation history obtained from the pulse- 
by-pulse processing had a significant amount of clutter. 
Consequently, a filtering-based technique had to be employed 


to extract the information component from the processed data. 
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Figure 3. LACE IQ Envelope, GMT DAY 200, 7603.875 
Seconds After Midnight 
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Figure 4. Power Spectrum of IQ Envelope, GMT DAY 200, 
7603.875 Seconds After Midnight 
This filtering-based technique consisted of taking the 
maximum value of a 2 second (125 point) moving histogram to 
determine the most likely Doppler separation track. Figure 5 
for 13 July 1990 and Figure 6 for 18 July 1990 illustrates the 


7 


results obtained by using this method. In analyzing Figures 5 
ana 6, it is evident that both histories are parabolic in 
nature. The biggest difference between the two Doppler 
separation histories is the existence of a flatter track for 
13 July. Obtaining the power spectrum of these Doppler 
separation histories will reveal the frequency components 


resulting from the structural dynamics of the craft. 
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Figure 5. Doppler Separation, GMT DAY 195, 2 Second Iteration 
(Fig 7a from Ref. 2] 
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Figure 6. Doppler Separation, GMT DAY 200, 2 Second Iteration 
(Fig 7b from Ref. 2] 


II. LACE SIGNAL PROCESSING 


LACE Signal Processing refers to the equations and 
procedures required to describe the signal processing of the 
received signal and the calculation of the power spectrum of 
this received signal. This process is represented in Figure 7, 
where LO represents the local oscillator and LPF represents 
the low pass filter. A/Dis the analog to digital converter. 
FFT represents the fast Fourier transform. S?(n) is the 


squared in-phase component and Sp (n) is the squared quadrature 


s,(n) 
o 
cos (w,¢) s?(n) 
y (k) 
PS 24s) 
sin(w,¢) s?(n) 
S 
So (n) 


Figure 7. LACE Signal Processing Block Diagram 
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component of the signal. The other terms in this figure are 
described in Equations (1) through (6). The equations used to 
describe the LACE Signal Processing were derived in References 
3 and 4. 

The signal received from the LACE satellite consisted of 
the components from both the forward boom's retroreflector and 
the body's retroreflector. Equation (1) is a mathematical 


representation of a noise free signal: 


S(t) =a,cos [(w,+W,,) E+,] +a,cos [ (W,+w,,) t+, ] (1) 


where the following terms apply to this equation: 
s(t)-is the received signal, 
a,-is the amplitude of the leading retroreflector, 
a,-is the amplitude of the body retroreflector, 
W,7is the carrier frequency, 


Wj, is the Doppler frequency of the leading 
retroreflector, 


Wyo-is the Doppler frequency from the body 
retroreflector, 


$,-is the phase return of the leading retroreflector, and 
g,-is the phase return of the body retroreflector. 
This signal is then processed by the laser radar to yield the 
in-phase and quadrature components of the signal. The in-phase 


and quadrature components are passed through a LPF and an A/D 


ia. 


converter. The two resulting signals are described 


mathematically by Equations (2) and (3): 


S,(n) = cos (w@g:NT+d,) +cos (wgnT+d,) (2) 


s_(n) = sin(o,,nT+d,) +sin(o,nT+d,) (3) 


where the following terms apply to these equations: 

S,(n)-is the in-phase component and 

S,(n)-is the quadrature component. 
The IQ envelope was obtained by adding together the squared 
in-phase and squared quadrature components generated by the 


radar. This signal is represented by Equation (4): 


2 2 
(ay+az)n i (a,a,)n 


$3 (n) = 4 2 


cos [ (@g,-g2) N+, -9,] Se) 


where the following term applies to this equation: 
si(n)-is the IQ envelope. 
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Equations (5) and (6) are used to calculate the power spectrum 


of the IQ envelope: 


Nea -4( 2% 
J (>) nk 


y(k) => s3(n)e (5) 
n=0 


PS53(k) =y (k) ¥(k) * 6) 


where the following terms applies to these equation: 


y(k)-is the complex signal (Fourier transform of IQ 
envelope), 


y(k)*-is the complex conjugate of the complex signal, 
N -is the transform length, and 
PS,2(k)-1s the power spectrum. 
This results in one peak located at the Doppler separation 


frequency as was indicated in Figure 4. 
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III. KALMAN FILTER 


The Kalman filter was developed in 1960 by Dr R. E. 
Kalman as an optimal recursive filter for the estimation of a 
state vector from measurement data corrupted by noise. It 
offered advantages over other filters such as the Wiener 
filter in that it reduces the mathematical complexity of the 
processing of large data strings. A block diagram of the 
Kalman filter is depicted in Figure 8, where DELAY represents 
one discrete-time delay. The rest of the terms in this figure 


are discussed in Equations (7) through (16). 


X(k+1|kK+1) 






X(k|k) 


Z(k+1|k) X(k+1|k) 


Figure 8. Kalman Filter Block Diagram 
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The Kalman filter equations are derived in References 5 
and 6. It is assumed that the process is time invariant. A 
mathematical model of the system is given below. The 


discrete-time state equation is represented by Equation (7): 


X(k+1) =FX(k) +Gv(k) (2) 


The measurement equation is represented by Equation (8): 


Z(k+1) =HX(k+1) +Dw(k) (8) 


where the following terms apply to the two preceding 
equations: 

X(k+1)-1s the updated state vector, 

X(k)-is the state vector, 

F -is the state transition matrix, 


v(k)-is the discrete time process noise assumed to be a 
zero-mean, white, random sequence, 


Z(k+1)-1s the measurement vector, 
H -is the gain through which the output leaves the system, 


w(k) -is the discrete time measurement noise assumed to be 
a zero-mean, white, random sequence, 


G -is the gain through which the process noise enters the 
system, and 
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D-is the gain through which the measurement noise enters 
the system. 


The equations involved specifically in the Kalman filter 
algorithm are discussed by dividing then into the prediction, 
innovation, gain and correction parts. In the prediction 
section, the conditional mean of the state vector (Equation 
9), the conditional state error covariance matrix (Equation 
10), and the predicted measurement vector (Equation 11) are 


computed: 


X(k+1|k) =FX(k|k) (9) 
P(k+1|k) =F P(k|k) F7™+GQG? (10) 
2(k+1|k) =HX(k+1|k) (15 


where the following terms apply to these three equations: 
X(k+1|k) -is the conditional state estimate, 
X(k|k) ~is the previous state estimate, 


P(k+1|k) -is the conditional predicted state error 
covariance matrix, 


Q -is the covariance of the process noise, 
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P(k|k) -is the previous predicted state error covariance 
matrix, and 


Z(k+1|k) -is the estimated measurement. 
The innovation section calculates the error between the 
measurement equations and the innovation covariance. Equation 
(12) calculates the measurement residual or innovation between 


Equations (8) and (11): 


e(k+1|k) =Z(k+1) -Z2(k+1|k) (12) 


Equation (13) determines the innovation covariance: 


S(k+1|k) =HP(k+1|k) H7+DRD? (13) 


where the following terms apply to these two equations: 
e(k+1i|/k)-is the innovation or measurement residual, 
S(k+1|k)-is the innovation covariance, and 
R-is the measurement noise covariance. 

The Kalman filer gain or weighting factor is found by Equation 


(14): 


K(k) =P(k+1|k) H™S(k+1|k) 7 (14) 
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Wheres 

K(k)-is the Kalman filter gain. 
The Kalman gain is the weighting factor that is placed on the 
measurement residual and the covariance prediction in the 
correction phase. Equation (15) corrects the old estimated 


state and the covariance correction is found by Equation (16): 


X(k+1|k+1) =X(k+1|k) +K(k) e(k+1|k) (15) 


P(k+1|k+1) =[I-K(k) H] P(k+1|k) (16) 


where the following terms apply to these two equations: 
X¥(k+1|k+1) -is the updated state estimate, 


fol cel) ceatl)) eeu the updated predicted state error 
covariance, and 


I -is the identity matrix. 
The preceding equations are then used recursively in the 
discussed order to obtain estimates of the state at a given k. 

There are two major factors that can affect the 
performance of the Kalman filter (Ref. 7], the first being the 
Kalman filter parameters such as process noise covariance, 
measurement noise covariance, and the initial conditions. 
These parameters are the fine tuning mechanisms of the filter. 
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It is seen in the Kalman filter equations that the gain is 
dependent on the prediction covariance and the measurement 
noise covariance. The prediction covariance is also dependent 
on the process noise covariance. If the process noise 
covariance in Equation (10) is increased, the prediction 
covariance increases. Therefore, the Kalman filter gain in 
Equation (14) can be considered as a trade off between the 
covariance of the process noise to the covariance of the 
measurement noise. With this in mind, as the process noise 
covariance increases, the Kalman filter gain increases and 
consequently, the bandwidth increases. This forces a faster 
transient response which leads to more noise in the estimates 
generated by Equation (11). By decreasing the measurement 
noise covariance in Equation (13), the same effect can be 
achieved. If the process noise covariance is decreased then 
the opposite effect will occur, which means that less noise 
will be present in the estimated states. With respect to the 
initial conditions chosen, the only part of the algorithm 
affected will be the transient part. As more data is processed 
the initial conditions fade eventually reaching a steady state 
value. By choosing a large prediction covariance more emphasis 
will be put on the measurements and less on the model in the 
transient phase. The second major factor affecting the Kalman 
filter performance is the model type. Since the model is used 
to generate the estimated states it should be as close to the 


physical phenomenon as possible. If the type of model chosen 


=o 


for a particular process is correct with only time constants 
slightly off, some degeneration will occur. On the other hand, 
if the system is modeled incorrectly, a model mismatch will 
result. A model mismatch can cause the estimate states to 
diverge from the actual states. It can be seen through the 
preceding discussion of the filter parameters and the choice 
of the model, the importance of correct modeling in the 
achievement of optimal performance from the Kalman filter 


(Ref. 8]. 
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IV. FIXED INTERVAL OPTIMAL SMOOTHER 


The Rauch-Tung-Striebel fixed interval optimal smoother 
was designed to be a post processing algorithm to be used in 
conjunction with a Kalman filter. This algorithm will improve 
the results obtained from the Kalman filter by utilizing the 
future information not available during the Kalman filtering 
process. The fixed interval optimal smoothing algorithm 
recalculates each estimate generated from the Kalman filter 
based on the information obtained for the entire set of data 
analyzed. This procedure generates what is called the smoothed 
estimates as seen in the block diagram of Figure 9, where the 
term ADV represents one discrete-time advance. The other 


terms are discussed in Equations (17) and (18). 





Figure 9. Fixed Interval Optimal Smoother Block Diagram 
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The Rauch-Tung-Striebel fixed interval optimal smoothing 
algorithm was derived in Reference 5. Equation (17) 


calculates a weighting factor: 


A(k) =P(k|k) F™P(K+1|k) Gi 


where: 
A(k)-is the smoothing algorithm gain, 


F -is the state transitional matrix from the Kalman 
filter, 


P(k|k) -is the prediction covariance from the Kalman 
filter, and 


P(K+1|k) -is the conditional prediction covariance from 
the Kalman filter. 


This weighting factor does not depend on past gains, just the 
conditional prediction error covariance and the previous 
prediction error covariance from the Kalman filter for a 
particular discrete-time. If more uncertainty exists in the 
forward filter, the weighting factor becomes larger. Equation 
(18) is the second equation involved in this algorithm which 


calculates the smoothed estimates: 


X(k|n) =X(k|kK) +A(k) (X(k+1|n) -X(k+1|k) ] (18) 
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where: 
n-is the final time, 
X(k|n) -is the smoothed state estimate, 
X(k|k) -is the state estimate from the Kalman filter, 
X(k+i|n) -is the previous smoothed state estimate, and 
X(k+1|k) -is the conditional state estimate from the 
Kalman filter. 
The next equation is not involved in the algorithm, but is 
useful in determining how well the smoothing is being 


accomplished. 


P(k|n) =P(k|k) +A(k) [P(k+1|n) -P(k+1|k)] A(k)? (19) 


where: 

P(k|n) -is the smoothed state error covariance 

P(kt+i|n) -is the previous smooth state error covariance 
These equations, excluding Equation (19), are used recursively 
in the discussed order to obtain the smoothed estimates. 

It is seen from the preceding equations that this 
algorithm has no performance parameter that can be adjusted. 
Consequently, it depends solely on the accuracy of the Kalman 
filter's parameters. The gain in Equation (17) is dependent on 
the covariance error of the previous and conditional values. 
This gain is then applied to the difference between the smooth 
estimate calculated for the previous data point and the 


23 


corresponding predicted value generated by the Kalman filter. 
The Kalman filter estimate is then adjusted by this weighted 
difference as in Equation (18) (Ref. 9]. This means that the 
predicted states, corrected states, predicted covariance, and 
corrected covariance for each data point must be saved during 


the forward processing operation. 
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V. KALMAN SMOOTHER APPLIED TO LACE 


The Doppler separation history for the LACE satellite is 
determined by utilizing the Kalman filter and the 
Rauch-Tung-Striebel fixed interval optimal smoother. Numerous 
preliminary steps were required to extract the Doppler 
separation from each pulse and ultimately obtain the Doppler 
separation history. The data first had to be read from the 
tape and converted into data files to be loaded into the 
Matlab environment for processing [Ref. 10}. An algorithm had 
to be developed to perform the preliminary processing of the 
data as discussed in the second chapter. The Kalman filter and 
the Rauch-Tung-Striebel fixed interval optimal smoother 
algorithms had to be implemented. A dynamic tracking window 
controlled by the Kalman filter had to be designed to track 
the desired frequency contained in the power spectrum for each 
radar pulse. Bad data had to be identified and eliminated 
during the processing phase. Adequate plotting routines had 
to be developed so that a good visual pulse-by-pulse analysis 
of selected data could be performed to ensure that the results 
were, in fact, what was expected. The following paragraphs 
discuss the previously-mentioned items in more detail with 
respect to the subroutines that were developed to perform 


these tasks. 
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The data received for 13 and 18 July 1990 consisted of 
approximately 9200 records per tape (120 seconds of data in 
binary form). Only 4500 records from each tape were processed 
in order to retain continuity with respect to the amount of 
data processed with the histogram filtering based method. 
Each record represents one radar pulse, which is 3.4 
milliseconds of IQ data digitized at 1.2 MHz. Subroutine 
DFSIOC found in the appendix, performs the task of extracting 
this information from the tape and coordinating the other 
subroutines to process the data. This subroutine enables radar 
pulses or records to be read from any specified run time to 
any other specified run time by adjusting the skip and count 
on the tape drive control line (!rsh srvl1 "dd if=/dev/nrmto 
ibs=16384 count=100 skip=0O > bin.dat). 

Subroutine BIN2INT found in the appendix was developed 
using Fortran [{Ref. 11] to convert the binary in-phase, 
quadrature, and timing information from the tape into an 
integer format capable of being loaded into the Matlab 
environment. Fortran was used in this application because 
Matlab does not have this conversion capability. The binary 
data is converted by BIN2INT, one record (16384 bytes or 
characters) at a time, into an integer format. Of the 16384 
characters, 1 to 16360 characters contained the IQ data for 
the radar pulse, 16361 to 16367 characters contained the 
timing information for the radar pulse, and the remaining 18 


characters are not used. The storage of the IQ data on the 
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tape alternated between a byte of in-phase and a byte of 
quadrature data. The timing information was stored in 
consecutive nibbles. The IQ information that subroutine 
BIN2INT converted is then stored in a data file called 
IQ DAT.dat and the timing information is stored ina data file 
called TimDat.dat. These data files are then sequentially 
loaded into the Matlab environment. 

Once the IQ data was converted into integer format and 
loaded into the Matlab environment, the next step was to 
process the data. The first, second, and third steps in 
subroutine DFSCAL, found in the appendix, performed the 
procedure discussed in the second chapter. These steps 
consisted of adding together the squared in-phase and squared 
quadrature components as in Equation (4) to produce the IQ 
envelope. The IQ envelope of a sample radar pulse is plotted 
in Figure 3. The conversion from the time domain to the 
frequency domain is accomplished by taking the fast Fourier 
transform of the IQ envelope and then calculating the power 
spectral density as in Equations (5) and (6). The results are 
shown in Figure 4 for the sample radar pulse depicted in 
Figure 3. 

The next step involved tracking the particular frequency 
of interest (the Doppler separation) contained in the power 
spectrum for each radar pulse. A dynamic tracking window 
controlled by a Kalman filter was designed to perform this 


function. In the development of the dynamic window, two 
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aspects had to be resolved: location and size. The location of 
the window was determined by centering it around the Doppler 
separation estimate obtained from the Kalman filter. The size 
of the window had to be large enough to track the change in 
frequency, but small enough to reject unwanted frequencies. 
Tracking these unwanted frequencies would cause the window's 
location to drift from the desired frequency. Window size 
control was achieved by setting the window's size equal to the 
prediction state error covariance from the Kalman filter plus 
two times the frequency resolution of the fast Fourier 
transform. The use of the prediction state error covariance 
plus two times the frequency resolution was empirically 
determined to be the parameter that produced the best results. 
The Doppler separation is then determined by obtaining the 
corresponding frequency of the peak magnitude within the 
window. This corresponding frequency of the peak magnitude is 
entered in Equation (12) of the Kalman filter as the 
measurement. The measurement is processed with the Kalman 
filter as discussed in the third chapter. 

The first obstacle encountered in processing the data was 
the clutter mentioned in References 1 and 2. This clutter or 
bad data can cause the tracking mechanism in subroutine DFSCAL 
to drift from the actual Doppler separation track if not 


addressed. Figure 10 was obtained by plotting the peak 
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Figure 10. Fixed Window, Unfiltered Data, GMT DAY 195 


magnitude within a fixed window from 4 kHz to 14 kHz. Each dot 
in this figure represents one radar return that could contain 
either the accurate Doppler separation information or the 
apparent clutter. In analyzing the power spectrums of selected 
radar pulses there exists three distinct types of radar 
returns: ones that contained accurate Doppler separation 
information; ones that contained accurate Doppler separation 
information corrupted by noise; and ones that are just noise 
or an extremely weak signal embedded in noise. The IQ envelope 
of a sample radar return that contains accurate Doppler 
separation information is depicted in Figure 11. Figure 12 
shows the power spectrum of this radar return and the dynamic 


tracking window. The dashed lines represent the window. The 
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solid lines indicate the power spectrum of the IQ envelope. 
The dotted line traces the outline of the power spectrum. The 
power spectral density being centered around a_ single 
frequency in Figure 12 is a good indication that this radar 
return contains accurate Doppler separation information. A 
sample radar return which contains accurate Doppler separation 
information corrupted by noise is represented by Figure 13 and 
its power spectrum is shown in Figure 14. In this case, the 
Kalman filter's filtering mechanism will eliminate any large 
deviations from occurring in the Doppler separation history. 
The third type of radar return is when there was just noise or 
an extremely weak signal embedded in noise. This type is 
represented in Figure 15. The power spectrum of this indicates 
the presence of noise as seen in Figure 16. Comparing Figure 
16 to Figure 14, a difference in magnitudes is revealed. This 
type of superfluous data is eliminated by placing a magnitude 
constraint on the particular frequency of interest within the 
dynamic window. The incorporation of the conditional statement 
(if MAGD(K+1)<MAGD(K)+1e2} MAGD(K+1)>MAGD(K)-1le2) into 
subroutine DFSCAL performs this task. This conditional 
statement differentiates between the magnitude of the last 
good radar return and the present radar return to determine if 
the relative magnitude is within a fixed distance of the last 
radar return. If this condition is met, the Kalman filter 


algorithm is used as described in the third chapter. On the 
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Figure 13. IQ Envelope Corrupted by Noise 
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Figure 14. Power Spectrum Corrupted by Noise 
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Figure 16. Amplified Noise in Frequency Domain 
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other hand, if this condition is not met, then the estimate of 
the Doppler separation frequency is not updated and the 
prediction covariance is increased by setting it equal to the 
conditional predicted state error covariance. This will cause 
the dynamic window to increase in size and place more emphasis 
on the next valid measurement within the window. 

The next obstacle encountered in the implementation of 
the Kalman filter, was determining the type of model best 
suited for the data. In developing the model for the Kalman 
filter, the histogram-based tracks for 13 and 18 July in the 
first chapter were analyzed to determined if a differential 
equation could be applied to the trajectory of the LACE 
satellite. After examining these figures, the development of 
an elaborate model to simulate the parabolic trajectory of the 
LACE satellite was virtually impossible to realize. The reason 
was due to the different elevations and azimuths the satellite 
could incur each time it was tracked. This leads to the 
slightly different trajectories as noted earlier with Figures 
5 and 6. A very simple model approach was then used. Adequate 
results were obtained using a first order model whose state is 
the Doppler separation in the frequency domain. Having the 
state of the model be a scalar quantity greatly simplifies the 
algorithm. To determine the value for the state transition 
matrix F, the parabolic nature of the histogram based tracks 
were taken into account. This meant that a frequency increase 


would occur at the start of a tracking run. After a maximum 
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value was reached during the tracking run, there would occur 
a frequency decrease. The only eigenvalue for the F that would 
not cause degeneration to occur, is the value one. If an 
eigenvalue greater or less than unity is used, the estimated 
frequency will diverge either for the increasing Doppler or 
the decreasing Doppler, depending on which value was chosen. 

When the preliminary processing of the data and the 
Kalman filtering had been accomplished, the Doppler separation 
estimates were ready to be smoothed. Subroutine DFSBFS found 
in the appendix was developed to perform this task. This 
subroutine implements the Rauch-Tung-Striebel fixed interval 
optimal smoother equations as discussed in the fourth chapter. 
Subroutine PLOT4 found in the appendix is utilized to plot the 
results obtained from this subroutine. 

The final step was to determine the value for the 
measurement noise covariance (R), the process noise covariance 
(QO), and the initial conditions that would optimize the 
processing algorithm. The initial conditions consisted of the 
prediction estimate covariance (P(0)), the Doppler separation 
estimate (XFREQ(0)), and the relative magnitude of the Doppler 
separation (MAGD(0)). Subroutine DFSVAL found in the appendix 
is used to initialize the parameters for all the subroutines. 
The initial Doppler separation frequency is determined by the 
point at which the analysis is started for a particular 
tracking run. The first radar return to be processed for 18 


July 1990 was depicted in Figure 3. The frequency of 12.012 
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kHz obtained from this figure is used to initialize the 
Doppler separation estimate for the filter. The prediction 
estimate covariance was set to 292.969 Hz, since this was the 
frequency resolution obtained. The initial magnitude is 
determined to be 100 by checking the beginning data. The 
value for the process noise covariance (Q) was varied from 
292.969 Hz to 7.324 Hz. The measurement noise covariance (R) 
was varied between 292.969 Hz and 2929.690 Hz. In empirically 
obtaining the optimal values for the process noise covariance 
and the measurement noise covariance, they were initially set 
at 292.969 Hz. By setting the process noise covariance and the 
measurement noise covariance to this value, the Kalman filter 
tracked every deviation no matter how off track they were. 
Figure 17 shows this undesirable result. The process noise 
covariance was when reduced to 146.485 Hz with the measurement 
noise covariance still set at 292.969 Hz in an attempt to gain 
more filtering from the Kalman filter. This showed some 
improvement in the noise rejection performance of the filter 
as indicated by Figure 18. Increasing the measurement noise 
ten times to 2929.960 Hz with the process noise covariance 
left at 146.485 Hz vastly improved the performance of the 
filter, which is seen in Figure 19. Based upon the previous 
results, the process noise covariance was further reduced to 
7.324 Hz leaving the measurement noise covariance at 2929.690 


Hz. A good track of the Doppler separation history was 
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Figure 19. Doppler Separation, GMT DAY 200, 
Q= 146.485, R = 2929.690 


obtained as indicated in Figure 20 for 18 July 1990. Figure 21 


shows the Doppler separation history for 13 July 1990. 
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Figure 20. Doppler Separation, GMT DAY 200, 
QOQ= 7.324, R = 2929.690 
20 


18 
16 
14 
12 


10 





1.1 1.20% 1.202 feos” loose Terooa ssi. 106 1.107 


GMT SECONDS AFTER MIDNIGHT x106 


Figure 21. Doppler Separation, GMT DAY 195, 
OQ = 7.324, R = 2929.690 
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VI. CONCLUSION 


The results obtained in the previous chapter have 
demonstrated the effectiveness of the Kalman filter in 
conjunction with the Rauch-Tung-Striebel fixed interval 
optimal smoother as a post processing utility in determining 
the Doppler separation history from the data. To compare the 
effectiveness of the two filtering-based techniques, Figure 22 
for 13 July 1990 and Figure 23 for 18 July 1990, were 
developed. Figure 22 is the combination of Figures 5 and 21, 


where Figure 23 is the combination of Figures 6 and 20. 
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Figure 22. Doppler Separation Composite, GMT DAY 195 
(Fig 7a from Ref. 2] 
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In both figures, the dashed lines are the results obtained 
from the histogram filtering-based method, the solid lines are 
for the method utilized in this thesis, and the dots represent 
the unfiltered data. Analysis of these figures demonstrates 
that the method used in this thesis of determining the Doppler 
separation history for both days is far superior to that of 


the histogram filtering based technique. 
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Figure 23. Doppler Separation Composite, GMT DAY 200 
{Fig 7b after REf. 2] 
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APPENDIX 
SUBROUTINES DEVELOPED TO PROCESS THE DATA 


The following subroutines are designed to function 
together with subroutine DFSIOC as the primary subroutine. The 
data is read from the tape in equally-sized blocks that are 
consecutively processed. The block sizes can vary from 1 to 
whatever size the computer is capable of handling. The number 
of blocks or groups can be varied from 1 to the amount of data 
to be processed. To use this package, the following must 
reflect the same values: the number of groups (GRNO) in 
subroutines DFSIOC and DFSVAL; the terminal count in the main 
loop of BIN2INT and the number of records per group (PSNO) in 
DFSVAL; the count in the tape drive control line of DFSIOC and 
the number of records per group (PSNO) in DFSVAL. These 
subroutines at present are set up to process 4500 records in 
groups of 100 (GRNO=45, PSNO=100). 


e DFSIOC 

To initiate this subroutine to process the data, just 
enter DFSIOC in the Matlab environment. DFSIOC will perform 
the task of extracting the information from the tape and 
coordinating subroutines DFSVAL, BIN2INT, DFSCAL, DFSBFS, 
PLOT3, and PLOT4. The data extracted from the tape is stored 
in a temp file called Bin.dat. It loads the data that was 
converted (binary to integer) from the temp file (Bin.dat) by 
BIN2INT into the Matlab environment for further processing. 
PSNO in this subroutine must reflect the same end count as 
PSNO in subroutine DFSVAL. In the tape drive control line, 
count is the number of records to be read and skip is the 
number of records to be passed over. Count must reflect the 
same value as PSNO in subroutine DFSVAL. The other options 
available are: to plot the data using PLOT3 and PLOT4; to save 
the data with the save function. 


t SUBROUTINE DFSIOC 
adfsval; SCALLS SUBROUTINE DFSVAL 
for GRNO=1:1:45 RMAIN LOOP TO PROCESS DATA 
{rsh srvl “dd ifz/dev/nrmt0O ibs=16384 &TAPE DRIVE CONTROL (READS) 
count=100 skip=0" >bin.dat 
{bin2int RCALLS SUBROUTINE BIN2ZINT 
firm /staff/thorngre/bin.dat R@CLEARS TEMPORARY DATA FILE 
load IQ _DAT.dat %*LOADS IQ DATA FILE 
load TimDAT.dat &% LOADS TIMING INFORMATION 
TimDAT1=[TimDAT1;TimDAT(:,6)]; SSAVES TIMING INFORMATION 
DAY=TimDAT(1,1); RSAVES GMT DAY OF TRACK 
dfscal; RS8CALLS SUBROUTINE DFSCAL 
clear IQ DAT TimDAT RCLEARS UNNECESSARY DATA 
end SEND OF THE MAIN LOOP 
dfsbfs RCALLS SUBROUTINE DFSBFS 
{rsh srvl mt rw RREWINDS TAPE TO BEGINING 
plots; %*CALLS SUBROUTINE PLOT3-OPT 
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plot4; 

save DSF.dat DSF -ascii 
*save XDSF.dat XMFREQ -ascil 
*save BXDSF.dat BXDSF -ascil 


tsave TimDAT1].dat TimDAT1] -aScil 


esave DAY.dat DAY -ascil 


e@ DFSVAL 


t*CALLS SUBROUTINE PLOT4-OPT 
@SAVES DATA IN ASCII FORMAT 


Declares and sets the parameters used in the subroutines 


DFSIOC, DFSCAL, PLOT1, 


and PLOT4. PSNO is the 


number of records in the group (GRNO). These must reflect the 
values used in DFSIOC and BIN2INT. Subroutine DFSVAL is called 


by subroutine DFSIOC. 


% 

GRNO=45; 

PSNO=100; 

RECLN=8160; 
Freg=292.969.*(0:2047); 
Time=(0:4079)/(1.2e6); 
TimDAT1=[]? 
G=zeros(PSNO*GRNO+1 ,1); 
P=zeros(PSNO*GRNO+1,1); 
PP=zeros(PSNO*GRNOt1,1); 
DSF=zeros(PSNO*GRNO+1,1); 
DSFW=zeros (PSNO*GRNOt+1,1); 
XDSF=zeros (PSNO*GRNO+1,1); 
XMFREQ=zeros(PSNO*GRNO+1 ,1); 
XFREQ=zeros (PSNO*GRNO+1 ,1); 
BXDSF=zeros (PSNO*GRNO+1 ,1); 
MAGD=zeros(PSNO*GRNO+1] ,1); 
DAY=(1,1); 

magfi=(1,1); 

fregfi=(1,1); 

magdi=(1,1); 

fregdi=(1,1); 
fregli=(1,1); 
fregui=(1,1); 

F=[{l]; 

G=(1]; 

D=[1]; 

C=(1l]; 

O0=[292.969*0.025]; 
R=[292.969*10]; 

P(1)=(292 .969*1]; 
MAGD(1)=[2.5e21]; 
XFREQ(1)=[292 .969*41]; 


@ BIN2INT 


SUBROUTINE DFSVAL 


tNO OF GROUPS OF PULSES 

%*NO OF PULSES PER GROUP 
t*LENGTH OF DATA PER PULSE 
RCREATES FREQUENCY VECTOR 
%*CREATES TIME VECTOR FOR IQ“%2 
SUNDETERMINED TIMIG VECTOR 
t*SETS UP KALMAN GAIN VECTOR 
*SETS UP COV PRED VECTOR 
RSETS UP CON COV PRED VECTOR 
SSETS UP DSF FIXED WINDOW VET 
RSETS UP DYNAMIC WINDOW VET 
*SETS UP ESTIMATED DSF VET 
@SETS UP CON EST STATE VET 
@SETS UP ESTIMATED STATE VET 
SETS UP SMOOTH EST DSF 
*SETS UP MAG VECTOR FOR D.W. 
%*GMT DAY OF TRACKING RUN 
t*MAX MAG IN FIXED WINDOW 
tINDEX OF MAX FREQ IN F.W. 
t*MAX MAG IN DYNAMIC WINDOW 
tINDEX OF MAX FREQ IN D.W. 
INDEX OF LOWER LIMIT D.W. 
t*INDEX OF UPPER LIMIT D.W. 
%TRANSITIONAL STATE MATRIX 
t*INPUT WEIGHTING FACTOR 
t*INPUT WEIGHTING FACTOR 
tOUTPUT WEIGHTING FACTOR 
PROCESS NOISE COVARIANCE 
t*MEASUREMENT NOISE COVARIANCE 
*INITIAL COV PRED ESTIMATE 
*INITIAL MAGNITUDE OF PSD 
INITIAL FREQUENCY ESTIMATE 


Has been written in Fortran to convert the binary data 
obtained from the tape into integer format. The input format 
is in binary byte form from the temp file Bin.dat. There are 


16384 bytes of data in a record of which 16360 bytes are 

inphase and quadrature data, 6 bytes are the timing 
information, and the rest are not used. There are no 
delimiters between values. The output format consists of I6 
format for the data and is stored in IQ DAT.dat. The output 
for the timing information is F4.4 format and is stored in 
- TimDAT.dat. The main loop in this subroutine must specify the 
numbers of records that are going to be processed per group 
(PSNO) in subroutine DFSVAL. BIN2INT at present is set up to 
process 100 records. This subroutine is called by subroutine 


DFSIOC. 


c SUBROUTINE BINZINT 
Cc DECLARES VARIABLES 
character*] cdat(16384) 
integer*2 jdat(8&192) 
integer*2 nibbles(0:11) 
integer*2 tcdat(0:2) 
real day,hour,min,sec,mil,TIME 
equivalence (cdat,jdat) 
S OPENS THE DATA FILE CALLED ‘’BIN.DAT’ FOR DIRECT, UNFORMATTED READ 
c AND OUTPUTS TWO FILES CALLED ‘IQ DAT.DAT’ AND ‘TIMDAT.DAT’. 
open(1,file=’/staff/thorngre/bin.dat’, access=‘direct’, 
& recl=16384,form=’unformatted’) 
open(2,f1le=’IQ DAT.dat’) 
open(3,file=’TimDAT.dat’) 
Cc LOOP TO READ DATA INFORMATION AND OUTPUT IT TO A FILE 
15 do 60, 1=1,100 
read(l, rec=1)cdat 
do 20, j=1,8192 


write(2,1020)jdat(j) 
1020 format (i6) 
20 continue 
c LOOP TO READ TIMING INFORMATION 


do 30, k=0,2 
tcdat(k)=jdat(8&181+k) 


30 continue 
e LOOP TO CALCULATE TIME FROM TIMING INFORMATION AND OUTPUT IT 
S TO A FILE 


do 40 1=0,11,4 
nibbles(l1+0)=and(rshift(tcdat(l1/4),12),15) 
nibbles(l+1)=and(rshift(tcdat(1/4),8),15) 
nibbles(l+2)=and(rshift(tcdat(1/4),4),15) 
nibbles (1+3 )=and(tcdat(l1/4),15) 

40 continue 
day=nibbles(0)*100+nibbles(1)*10+nibbles (2) 
hour=nibbles(3)*10+nibbles (4) 
min=nibbles(5)*l0+nibbles(6) 
sec=nibbles(7)*l0+nibbles(8) 
mil=nibbles(9)*l00+nibbles(10)*l10+nibbles(11) 
time=60*60*hour+60*min+sec+(mil/1000) 
write(3,1030)day,hour,min,sec,mil,time 

1030 format (£4,1x,£4,1x,£4,1x,£4,1x,£4,1x,£10.4) 

60 continue 

end 
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@ DFSCAL 

Processes records that contain the inphase and quadrature 
data to obtain the IQ envelope, power spectrum of the IQ 
envelope, and the estimated Doppler separation history. The 
Kalman Filter algorithm is implemented to perform two 
functions. The first is to provide control for the dynamic 
window's location and size. The second is to provide a good 
estimation of the Doppler separation from the measured data. 
The conditional statement is used to check the magnitude of 
the data. The fixed window in this subroutine provides an 
observation of the unfiltered data. The options available in 
DFSCAL are to view IQ envelope and the power spectrum of the 
IQ envelope subroutine PLOT1l or subroutine PLOT2. This 
subroutine is called by DFSIOC. 


% SUBROUTINE DFSCAL 
for K=(GRNO*PSNO-PSNO+1 ):1:(GRNO*PSNO) 
IQS NO=1; 
for k=1:2:RECLN 
IQS_DAT(IQS_NO)=(IQ_DAT(k+(IQS_NO-1) ... 
*8192))*2+(IQ_ DAT(k+ ... 


SLOOP TO CALCULATE DSF/PULSE 
NINITIAL COUNTER FOR IQ DAT 
*LOOP TO CALCULATE IQ*2 


% 
% 


1+(IQS_NO-1)*8192))*2; 
IQS_NO=IQS NO+1; 
end 
IQS_FFT=f£f£t(IQS_DAT,4096); 
IQS _PSD=IQS_FFT.*conj(IQS FFT); 
[magfi freqfiJ=max(IQS_PSD(15:49)); 
freqfi=fregfit+l4; 
DSF(K)=Freq(freqfi); 
XMFREQ (K+1)=F*XFREQ(K); 
XDSF (K)=C*XMFREQ (K+1); 
PP(K)=(F*P(K)*F’)+(G1*Q*G1’); 
G(K)=PP(K)*C’*inv(C* (PP (K)*C’+D*R*D’)); 
fregli=round( (XDSF(K)-PP(K))/292.969); 
frequi=round( ((XDSF(K)+PP(K))/292.969)+2); 
[magdi freqdiJ=max(IQS PSD(freqli... 
:frequi)); 
freqdi=fregli+freqdi-1; 
DSFW (K)=Freq(freqdi); 
MAGD (K+1)=magdi; 
if MAGD(K+1)<MAGD(K)+1le2 |... 
MAGD (K+1)>MAGD(K)-1e2, 
XFREQ (K+1)=XMFREQ (K+1)+G(K) 
* (DSFW(K)-XDSF(K)); 
P(K+1)=(eye(1)-G(K)*C)*PP(K); 
else 
XFREQ (K+1)=XMFREQ(K+1); 
P(K+1)=PP(K); 
DSFW(K)=0; 
end 
plotl 
plot2 


end 
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REND OF LOOP TO CAL IQ*2 
RCALCULATES FFT OF IQ*2 
%CALCULATES POWER SPECTRUM 
SFIXED WINDOW 4KHz AT 14KHz 
RS INDEX CORRECTION 

SDETERMINE DSF FROM INDEX 
SOPDATES ESTIMATED STATE 
SCALCULATE ESTIMATED DSF 
SUPDATE CON PRED COVARIANCE 
%CALCULATE KALMAN FILTER GAIN 
SLOWER INDEX FOR WINDOW 
SUPPER INDEX FOR WINDOW 
SINDEX MAX MAG WITHIN WINDOW 


% INDEX CORRECTION 
%DETERMINES FREQ FROM INDEX 
&SAVES MAG FOR EACH LOOP 
SCONDITIONAL TEST FOR PULSE 
%PULSE CONTAINS USEFUL DATA 
SUPDATES STATE 


SUPDATES COVARIANCE PREDICTION 
%PULSE HAS NO USEFUL DATA 
SSTATE EQUALS CON EST STATE 
S%SET COV PRED TO CON COV PRED 
NO MEASUREMENT OBTAINED 

REND OF CONDITIONAL STATEMENT 
%*CALLS SUBROUTINE PLOT1-OPT 
%CALLS SUBROUTINE PLOT2-OPT 
tEND OF LOOP TO CALCULATE DSF 


e DFSBFS 
Was developed based on the Rauch-Tung-Striebel fixed- 


interval optimal smoother algorithm to provide a smoothed 
estimate Doppler separation history. This subroutine is called 


by DFSIOC. 


% SUBROUTINE DFSBFS 
BXDSF (1:PSNO*GRNO+1 ) =XFREQ(1:PSNO*GRNO+1); SINITIALIZES 1ST VALUE BDSF 
for n=PSNO*GRNO:-1:1; &8LOOP TO CALCULATE BACKWARDS 

A=P (n+1)*F’*(inv(PP(n)))7 RCALCULATES SMOOTHING GAIN 
BXDSF (n)=XFREQ(n)+(A* (BXDSF(n+1) ... RUPDATS SMOOTH ESTIMATE DSF 
-XMFREQ(n+1)))7 
end SEND OF LOOP CAL BACKWARDS 
e PLOT1 


Plots the IQ envelope for the radar pulse to the monitor. 
The option in PLOT1 is to save the plot by using the meta 
function. This subroutine is called by DFSCAL. 


% SUBROUTINE PLOT] 
titlel=([’SECONDS AFTER MIDNIGHT ‘’,... STIMING INFORMATION HEADER 
num2str(TimDAT1(K,1))])3 

title2=([{'’GMT DAY ‘’,num2str(DAY)]); RTIMING INFORMATION HEADER 
plot(Time(1:4080),IQS_DAT(1:4080)),grid SPLOTS IQ ENVELOPE DATA 
title({’LACE IQ ENVELOPE,’,title2]) SPLOT TITLE 
text(1,1,titiel, “3c j; RPLOTS TIMING INFORMATION 
ylabel( ‘RELATIVE MAGNITUDE’), SPLOTS Y-AXIS LABEL 
xlabel(’TIME (ms)‘), %PLOTS X-AXIS LABEL 

meta sxsl SSAVES PLOT-OPTIONAL 

®@ PLOT2 


Plots the power spectrum of the IQ envelope and the 
aynamic window for the radar pulse to the monitor. The option 
in PLOT2 is to save the plot by using the meta function. This 
subroutine is called by DFSCAL. 


% SUBROUTINE PLOT2 
titlel=([{’SECONDS AFTER MIDNIGHT ’,... &STIMING INFORMATION HEADER 
num2str(TimDAT1(K,1))]); 
title2=([’GMT DAY ’,num2str(DAY)]); %TIMING INFORMATION HEADER 
plot(Freg(freqli-15:frequit15),IQS_PSD ... «PLOTS OUTLINE PWR SPECTRUM 
(fregli-15:frequit+l15),’:r’),hold on SRETAINS PLOT IN MEMORY 
YIGRA=O:magi:magi;Y11GRA=[magi magi]; SSETS UP DYNAMIC WINDOW VET 


X1GRA=Freg(freqli)*(ones(l:length(Y1IGRA))); 

X2GRA=(Freq(frequi)* (ones(l:length(YIGRA))))‘; 

Y2GRA=Freg(freqli):Freq(frequi)-Freg ... 
(freqli):Freg(frequi); 


plot (X1GRA, YIGRA, ‘--g',X2GRA,YIGRA, ’--g'), *PLOTS SIDES DYNANIC WINDOW 
plot(Y2GRA,Y1IGRA,’--g'), SPLOTS TOP OF DYNAMIC WINDOW 
magdis=[zeros(freqli-15:frequit15);... SSETS UP POWER SPECTRUM 
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IQS_PSD(freqli-I5:frequitl5)... 
pzeros(freqli-15:frequit+l5)]; 

indexs=(Freq(freqli-l15:frequi+15);Freq ... 
(freqli-l15:frequitl5);Freq(freqli ... 
-15:frequit+15)]; 

magdis=(magdis(:)']; 

indexs=[indexs(:)’]; 


plot(indexs,magdis,’-r’),grid %*PLOTS POWER SPECTRUM 
title([{’POWER SPECTRUM OF IQ ENVELOPE,’,title2]) 
bext(l1,1,titlel,’sc’); RPLOTS SECONDS AFTER MIDNIGHT 
xlabel(’FREQUENCY (KHz)’), RS@PLOTS X-AXIS LABEL 
ylabel(’RELATIVE MAGNITUDE’), SPLOTS Y-AXIS LABEL 

&meta sxs2 RSAVES PLOT-OPTIONAL 

hold off QRELEASES PLOT FROM MEMORY 

e PLOT3 


Plots the unfiltered data in the fixed window and the 
estimated Doppler separation history to the monitor. The 
option in PLOT3 is to save the plot by using the meta 
function. This subroutine is called by DFSIOC. 


% SUBROUTINE PLOT3 
title2([’GMT DAY ‘',num2str(DAY)]); STIMING INFORMATION HEADER 
title3(('Q = ',num2str(Q)]); SPROCESS NOISE COVARIANCE 
title¢({’R = ‘,num2str(R)]); %MEASUREMENT NOISE COV 


axis({TimDAT1(1,1) TimDATI1(length(TimDAT1)... %*DETERMINES AXIS FOR PLOT 
7,1) O 20000]) 

plot (TimDAT1,DSF(1:length(TimDAT1),1),’.r’,...%PLOTS DATA 
TimDAT!,XDSF(l:length(TimDAT1),1),’-g'); 


title([’DOPPLER SEPARATION,’,title2]) SPLOTS TITLE 
Eext(2,1,title3,’sc’), %PLOTS PROCESS NOISE Cov 
text(3,1,title4,’sc’), R%PLOTS MEASUREMENT NOISE COV 
xlabel(’SECONDS AFTER MIDNIGHT’), RQPLOTS X-AXIS LABEL 
ylabel(’DOPPLER SEPARATION (KHz)’), SPLOTS Y-AXIS LABEL 

% meta sxs4 RSAVES PLOT-OPTIONAL 

@ PLOT4 


Plots the smoothed Doppler separation history to the 
monitor. The option in PLOT4 is to save the plot by using the 
meta function. This subroutine is called by DFSIOC. 


% SUBROUTINE PLOT4 
title2({’GMT DAY ‘,num2str(DAY)]); STIMING INFORMATION HEADER 
title3([’Q = ’,num2str(Q)]); &PROCESS NOISE COVARIANCE 
title4([’R = ',num2str(R)]); &MEASUREMENT NOISE COV 


axis({TimDAT1(1,1) TimDAT1(length(TimDAT1)... %*DETERMINES AXIS FOR PLOT 
,1) 0 20000)) 


plot (TimDAT1,BXDSF(l:length(TimDAT1),1),... RPLOTS SMOOTH EST DSF 
Ja 4) yo Pablo Fe 
title([{ ‘DOPPLER SEPARATION, ‘’,title2]), QSPLOTS TITLE 
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text (2,0, tities, Se.), ®8PLOTS PROCESS NOISE CoV 


text(j,1,title4,’sce’), *PLOTS MEASUREMENT NOISE COV 
xlabel(’SECONDS AFTER MIDNIGHT’), RPLOTS X-AXIS LABEL 
ylabel(’DOPPLER SEPARATION (KHz)’), RPLOTS Y-AXIS LABEL 
% meta sxs¢4 SSAVES PLOT-OPTIONAL 
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